function D_new=osl_channelrepairMEG(S)

global OSLDIR

% GW 2013

D=S.D;

fname=fnamedat(D);

S2=[];
S2.D=D;
S2.newname=[path(D) filesep 'C' fname(1:end-4) '.mat'];
D_new=spm_eeg_copy(S2);

% three sensor types: grad orient 1, grad orient 2, mags
for iSensType = 1:3
    
    % prepare a layout, neighbourhood structure and connectivity matrix for
    % this subset of channels
    if iSensType == 1
        lay_fn = fullfile(OSLDIR,'layouts','neuromag306planar3.lay');
    elseif iSensType == 2
        lay_fn = fullfile(OSLDIR,'layouts','neuromag306planar2.lay');
    else
        lay_fn = fullfile(OSLDIR,'layouts','neuromag306mag.lay');
    end
    cfg=[];
    cfg.layout = lay_fn;
    lay = ft_prepare_layout(cfg);
    
    % get the indices into the channels in the D-object that contain this
    % sensor type
    Dlabels = D.chanlabels;
    [labels,inds_into_D] = intersect(Dlabels,lay.label);
    
    % neighbours
    cfg=[];
    cfg.method            = 'distance';
    cfg.neighbourdist     = 0.4; % quite big.  NB average is inversely weighted by distance.
    cfg.layout            = lay;
    cfg.channel           = labels;
    cfg.feedback          = 'no';
    neighbours            = ft_prepare_neighbours(cfg);
    
    % connectivity matrix
    cfg = [];
    cfg.neighbours = neighbours;
    cfg.channel = labels;
    connectivityMatrix  = channelconnectivity(cfg);

    bch_list=badchannels(D);
    [bch_list,bch_inds_reduced,~] = intersect(inds_into_D,bch_list);
    [gch_list,gch_inds_reduced]   = setdiff(inds_into_D,bch_list);
    A=chanlabels(D);
    chan_names_bad  = A(bch_list);
    
    connectivityMatrixGd = connectivityMatrix(:,gch_inds_reduced);
    if any(~(logical(sum(connectivityMatrixGd,2))))
        error('No good nearby channels!');
    end
    
    for iRepair = 1:numel(chan_names_bad)
        thischanindx_reduced = bch_inds_reduced(iRepair);
        thischanindx_intoD   = bch_list(iRepair);
        
        disp(['Repairing ' chan_names_bad{iRepair}])
        
        connected = connectivityMatrixGd(thischanindx_reduced,:);
        connected_inds_into_D = gch_list(find(connected));
        
        distance = sqrt(sum((lay.pos(gch_inds_reduced,:) - repmat(lay.pos(thischanindx_reduced, :), numel(gch_inds_reduced), 1)).^2, 2));
        recip_distance = (1./distance);
        recip_distance = recip_distance(connected)';
        recip_distance = recip_distance ./ sum(recip_distance);
        
        for iTrial = 1:D.ntrials
            dat = D(connected_inds_into_D,:,iTrial);
            newdat = recip_distance * dat;
            D_new(thischanindx_intoD,:,iTrial) = newdat;
        end
    end
    
end % for iSensType = 1:3

D_new.save;